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g '. Abstract 

Properties of weighted averages are studied for the general case that the individual 
measurements are subject to hidden correlations and have asymmetric statistical as 
well as systematic errors. Explicit expressions are derived for an unbiased average 
with a well defined error estimate. 
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§ 1 Introduction 

Combining independent measurements of a physical quantity in order to summarize the in- 
formation from different sources into one number is usually done by calculating a weighted 
>^ | average, where each measurement contributes with a weight which is inversely proportional to 

V ' its variance. 
Oh' 

In practical applications, however, one is often confronted with the problem that the mea- 
surement errors which are quoted define a 68% confidence level interval around a Maximum 
Likelihood estimate, whereas one would need an unbiased estimator with its standard devia- 
^ ■ tion in order to form a weighted average. In addition, one usually has to deal with statistical 
and systematic uncertainties, and different results may be correlated, e.g. due to common 
uncertainties in some theoretical prediction, while a quantitative estimate for the size of the 
correlation does not exist. 

This paper presents a systematic approach to form a meaningful weighted average under such 
conditions. After a short reminder of the properties of weighted averages, it first is shown how 
to deal with hidden correlations in the data. Then a detailed analysis of parameter estimates 
based on the Maximum Likelihood Method is performed, and the consequences when forming 
a weighted average are studied in detail. The results are illustrated by means of a numerical 
simulation. Finally, the issue of systematic errors and how to combine statistical and systematic 
errors is addressed. 



2 Weighted Averages 

Let {x^ be a set of n unbiased measurements of a parameter /i, with expectation values (x^ 
and covariance matrix Cij(x) given by 

(x^ = fi and Cij(x) = ((xi - fi)(xj - fi)) = (xiXj) - fi 2 . (1) 
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From {xi} a new estimator x for /i can be constructed, which has a smaller variance than any 
of the Xi, by minimizing 

V Y^i-r, ■ x)C^\x). (2) 

Note that the x 2 -function is constructed in such a way, that it is invariant with respect to 
arbitrary linear transformations of the variables with a regular transformation matrix. In case 
the individual Xi are independent, the covariance matrix becomes diagonal, Cij(x) = cxfSij and 
eq.(|]) simplifies to 

X =L ~2 • ( 3 ) 

The value x which minimizes x 2 an d thus in the least squares sense is optimally consistent with 
all measurements Xi is found to be 

x = ^WiXi with Wi = and ^ = E ~2 ' (4) 

and the variance of x is obtained by standard error propagation as 

C(x) = a\x) = WiWjdj = £ wlal = N' 1 . (5) 

ij i 

Another relevant quantity is the expectation value of x 2 at the minimum 

(xLn) = E ( ^2^ ) = E (J $- - N (x 2 ) = n - 1 (6) 

Note that although the value of (x 2 ) has been derived for the case of a diagonal covariance 
matrix, the result is valid in general since one always can find a linear transformation of variables 
which diagonalizes Cy while x 2 by construction is invariant. 

Minimizing J2i w2(J i with respect to subject to the constraint J2i w i — 1> one sees that the 
weights Wi used for the weighted average eq.@ lead to the smallest possible variance for x if 
the individual measurements are uncorrelated. 

In case that correlations are present, one can still calculate a weighted average according to 
eq.(|]), which always is an unbiased estimator for fi, even if it no longer is the value which 
minimizes eq.([|). It also does not have the smallest possible variance, although the loss in 
precision usually is only marginal. In practical applications one thus might prefer the robust 
average eq.(^) to the formally optimal result from minimization of eq.(^). As this requires 
to calculate the inverse of Cij, it can only be performed reliably if the off-diagonal elements 
of are known precisely. Otherwise numerical instabilities are likely to render the formally 
optimal result x meaningless. Note, however, that while it is perfectly meaningful to calculate 
a weighted average ignoring off-diagonal elements of the covariance matrix, the full matrix 
must be taken into account for the variance of the average. 

For uncorrelated inputs the weighted average eq.(Q) is the optimal way to combine measure- 
ments in order to minimize the variance of the result. It is independent of the actual shape of 
the probability density functions (PDFs) involved. All one needs to know are the values Xi and 
the variances of of the individual measurements. Note that the confidence levels associated 
to the variances may vary with the shape of the PDF. The weighted average, on the other 
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hand, has a PDF which is a convolution of the PDFs of the individual X{, and according to 
the central limit theorem with an increasing number of contributions quickly converges towards 
a Gaussian, where the range x ± u{x) corresponds approximately to an 68% confidence level 
interval. 

3 Hidden Correlations 

In practical applications one sometimes is confronted with a situation where measurements 
have to be averaged which are suspected to be correlated, but where the correlations are very 
hard to quantify. While, as argued before, a weighted average can still be calculated, it is 
not obvious how to give a meaningful estimate for the error of the estimate. One possible 
approach to this problem has been discussed in [p]]. The basic idea was to use the value of 
the x 2 calculated according to eq.(^) as an indicator and adjust a global correlation coefficient 
between all measurements such that it becomes equal to its expectation value. Here a similar 
argument will be presented, which is more convenient as it leads to a much simpler scaling 
prescription. 

A value \ 2 larger than its expectation value can result if anti-correlations are present, or if the 
measurement errors are underestimated. In case of the latter, iV -1 from eq.(|^) underestimates 
the variance of x, in the former case it is overestimated. In order to be conservative, a common 
practice |J is to assume that a large \ 2 is a consequence of too small measurement errors 
and apply a common scaling factor x 2 jin — 1) to all variances to obtain a x 2 identical to its 
expectation value. Although one certainly can question this approach, it has to be pointed 
out, that unless there are severe discrepancies in the data, this procedure still gives an error 
estimate for x which is smaller than any of the individual errors &i of the data. 

The other indication for a possible problem is a too small value for \ 2 ■ This can either mean 
overestimated measurement errors or positive correlations between the data. While the former 
would imply that the error of x in reality is smaller than iV -1 , the latter would require a larger 
error. Being conservative, again an estimate for the larger error shall be derived. As these 
considerations only affect error estimates, the calculations can be simplified without loss of 
generality by assuming (xj) = (x) = 0. This implies (x 2 ) = of and (x 2 ) = C(x), and the 
expectation value (xLin) ec i-@ simplifies to 

(xLn) = n-NC(x). (7) 

If the individual measurements are uncorrelated, then one has C(x) = iV -1 , and as shown 
before (x m in) — n — 1. If there are reasons to suspect that the measurements Xi are not 
independent, then one can re-interpret eq.(0), assuming that the actual best fit x 2_vame is 
equal to its expectation value and derive a variance estimate 

C(x) = ±(n- X 2 ). (8) 

This expression constitutes a simple scaling prescription for the case that the \ 2 °f a weighted 
average is smaller than expected. Like the scaling procedure for the errors in case of too large 
X 2 values, it defines an effective way to take imperfections or incomplete knowledge about the 
input data into account. 
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Like the weighted average, also the simple scaling eq.@ of the uncorrelated variance estimate 
N^ 1 is readily generalized to m-dimensional vector- valued measurements Xj. It will also be 
demonstrated explicitly that eq.([8]) is equivalent to the assumption of a fixed common covariance 
between all measurements. In contrast to this, the more complicated procedure discussed in 
jfj was based on the assumption of a common correlation coefficient. In practice both scaling 
schemes give rather similar results. Numerically the method |l] yields slightly smaller errors, 
indicating that the assumption of a common correlation coefficient tends to assign a higher 
weight for more precise measurements than the assumption of a common covariance. 

Ignoring correlations between measurements, the generalization of eq.(|3]) to vector valued mea- 
surements Xi with covariance matrices Cj is given by 



X 
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xfa 



> Xi 



X 



(9) 



Note that in this expression the index i enumerates entire vectors rather than individual compo- 
nents. Those vectors Xi are assumed to be uncorrelated; correlations between the components 
of Xi are described by the covariance matrices Ci. The minimum of eq.(O) is obtained for 



x 



N- 1 C^Xi with N = J2 Ci 1 - 



(10) 



Note that now x is a vector of dimension m. Generalizing eq. (|(J) the value of \ 2 at the minimum 
then can be expressed as 



xL^Tr^r 1 



Xixf — Nxx T 



(11) 



Without loss of generality one can again make the simplifying assumption (xj) = (x) = 0. This 
leads to (x^xf^ = Ci and (xx T ^ = C(x), and the expectation value of x min becomes 

= n ■ m - Tr (n (xx t )) = n ■ m - Tr (NC(x)) . (12) 

At this point one could directly apply the same short argument as before and derive the 
generalization of eq.@ for m-dimensional measurements. However, in order to illuminate a 
little bit the background of the scaling rule, here a different route shall be taken, making the 
explicit assumption that the correlation between the vectors X{ is such, that the covariance 
matrix Ci for each vector Xi is the sum of a common contribution C and a specific component 

Si 

j C for i ^ j 

a 



d = C + Si i.e. 



/ -> ->T 

{ XiXj 



for i = j 



(13) 



One then obtains 



C(x) = (xx T ) = N~ l 1 + NC - \ J2 C^CCr^j N' 1 



(14) 



where 1 denotes the unit matrix, and 



m ■ (n 



Tr 



NC 



Y,C^CC7An 



-1 



(15) 



For C = the covariance matrix of the average is given by iV 1 and one has (x m m) = m ' 
(n — 1). A significant deviation of the Xmm fr° m this expectation value can be taken as an 
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indication for correlations C between the individual measurements. The simplest way to take 
these correlations into account for the error of the average, is by scaling the uncorrelated 
estimate N~ l by a factor (1 + c). From eq.([T3|) then follows that C satisfies the relation 



NC 



Y t CZ x CCZ^N- 1 = cl. (16) 



For a given c this is a linear equation which can be solved numerically for C by iterating for 
example the fixed point condition 



C n+1 = N- 1 + N- 1 (j2 C^C n cA 



N- L (17) 



and setting C = c Coo. In practical applications this only has to be done if one wants to 
extract C in order to get a quantitative estimate for the size of the correlations between the 
measurements. Otherwise it is sufficient to fix the value of the parameter c by requiring the 
observed x 2_va hie to be equal to the expectation value (Xmm)- With eq. (|l6"D and eq. ([l5|) one 
obtains 

Xmin = m ( n - 1) - m C (18) 




C(x) = N- 1 [n-^\. (19) 



and a scaled variance estimate 



For m = 1 the result eq.(j8|) is recovered. 

To summarize, unless correlations between measurements are really well understood, it appears 
advisable to employ the weighted average ignoring all correlations as a robust procedure to 
combine several inputs into one number. Care, however, has to be exercised to give a realistic 
error estimate. In the absence of other information the x 2 °f the average is a useful indicator 
for potential problems. The conservative approach would be to interpret a x 2 smaller than its 
expectation value as evidence for positive correlations between the measurements, and a x 2 
larger than its expectation value as evidence that the errors are underestimated. In both cases 
the simple estimate N' 1 for the variance of the weighted average would have be to increased, 
which is conveniently done according to the following scheme 

for x 2 — m{n — 1) 
C/(.r): .V" l { V \ (20) 

for x 2 > rn(n — 1) 

It has to be emphasized, that such a scaling should only be performed if one has reasons to 
believe that the variance estimate iV -1 of the weighted average underestimates the true un- 
certainties, or if one wants to quote a conservative error. Otherwise, since the ^-distribution 
especially for a small number of degrees of freedom has a large relative width, one will system- 
atically bias the error to large values. 



4 Statistical Errors 

In order to form a weighted average of independent data one needs unbiased measurements 
with known variances. On the other hand, in order to quantify the statistical precision of a 
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measurement, usually a 68% confidence level interval is quoted, which is constructed such that 
the true value of the parameter is inside this interval in 68% of all cases. For a Gaussian 
distribution the 68% confidence level interval is given by x ± a, where x is the mean value and 
a the rms-width of the distribution, i.e. the square of the error is the variance. In general one 
needs to know the shape of the PDF in order to extract the variance from the error interval. 

If the shape of the PDF is not explictly specified, one has to resort to certain reasonable 
assumptions in order to proceed. For the following it will be assumed that in absence of other 
information the primary result of a measurement has a Gaussian PDF. Using the Maximum 
Likelihood Method this measurement then is interpreted in terms of a physical parameter, with 
a nominal result and an error range covering a 68% confidence level interval. Under these 
conditions a rigorous interpretation can be attached to asymmetric errors even if no further 
information about the PDF is given. 

4.1 Interpretation of Asymmetric Errors 

As explained above, a measurement x is assumed to scatter around a mean value (x) according 
to a Gaussian PDF with variance a 2 . The precision of the measurement a is assumed to be 
known, e.g. from first principles or some calibration procedure. The mean value (x) is related 
to a physical parameter \i via a function s(/x), (x) = s(/i). The likelihood function p{x\fi) to 
observe a value x given /z, then is 

In the framework of the Maximum Likelihood Method the function eq.(^TJ) is used to extract 
an estimate fi m i for \x from a given measurement x as the parameter which maximizes the 
likelihood. The solution is evidently given by [i m i = s~ 1 (x). Note that the likelihood function 
here is merely a tool to extract an estimate [i m i for \i by looking for the value which maximizes 
p(x\fj), i.e. no attempt is made to interpret the likelihood function as a PDF for fi. Note that 
in a Bayesian approach using a flat prior distribution, one could actually consider eq.(pT|) as a 
probability density for fi and take e.g. the average (/i) of p(x\p) instead of estimate 
for fi, although, as shown below, the Maximum Likelihood estimate /i m / usually is the better 
choice. 

In addition to the parameter estimate also an error interval has to be given. This is done by 
translating the error interval x + ] of the measurement to an error interval for the Maximum 
Likelihood estimate according to \p-,[i+] = [s _1 (a;_), By conservation of probability 

both intervals have the same probability content. If p(x) denotes the PDF for the measurement 
x, then the PDF for the Maximum Likelihood estimate q(/J* m i) is given by 

/oo 
dx p(x) 8{ii m i - s~ 1 {x)) (22) 
-oo 

and one obtains for the probability content 

/ dfi mi q(nmi) = dfi m i / dx p{x) S^Umi - s (x)) = dx p(x). (23) 

J fl— Js~ 1 (x-) J— OO JX — 

Usually the error interval for the measurement is given by [x-, x + ] = [x — a, x + a] which in case 
of a Gaussian PDF includes the expectation value (x) in approximately 68% of all cases. The 
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same confidence level then also applies for [//_,//+]. From eq.(|21|) one sees that this definition 
of error interval corresponds to finding regions where the logarithm of the likelihood function 
stays within half a unit of the maximum, i.e. it is the shortest interval which can be constructed 
for a given probability content. This concludes the discussion of what will be assumed to be 
the connection between a measurement x ± a and the corresponding Maximum Likelihood 
parameter estimate fx m i ± a±. 



It is now possible to study the expectation value and variance of Li m i- With eq.([22|) the expec- 
tation value becomes 

(fimi) = J dfimi fi m i qijimi) = J dx p(x) s~ 1 (x). (24) 

It is clearly not possible to determine this expectation value for arbitrary functions s~ 1 (x), 
but one can address the problem in a systematic way by expanding s _1 (x) around (x) = s(/x). 
Since an asymmetric error will not come about for linear functions s(fx), the expansion should 
at least go to second order. Assuming that for practical applications the non-linearities will be 
small, the expansion will be truncated there. 

s^ 1 (x) = fi + ay + (3y 2 with y = x — (x) . (25) 

Substituting this into eq.(|24]) and using (y 2 ) = a 2 then yields 

(fimi) =/i + /V. (26) 

One finds that the maximum likelihood estimate fx m \ is a biased estimator, with a constant 
bias proportional to the curvature of s~ l (x). An unbiased estimator can be constructed by 
exploiting the information contained in the asymmetric errors. As discussed before, the error 
range /i± = fi m i ± cr ± is defined by the condition fx± = s _1 (x ± er), which yields 

o+ = ao + 2f3oy + f3a 2 and o\_ = aa + 2(3ay - f3o 2 . (27) 



The difference of the two is proportional to the bias term in eq.(PB"D so that an unbiased estimator 
pL for fi is given by 

A = Pmi - ~ a -)- ( 28 ) 

At first glance this result may be surprising, since it implies that in case the positive error is 
larger than the negative one, the unbiased estimate (x for the parameter is in fact smaller than 
fi m i. This underlines the fact that the asymmetric error does not describe a likelihood function 
for /i. Instead, a large positive error means that the value Li m i on average will overestimate the 
true value, which is then compensated by subtracting half the difference of the errors. 

In terms of the expansion eq. (P5[) the unbiased estimate fx is given by 

fl = Li + ay + f3y 2 - f3a 2 (29) 

and with 

fe " )= V^/^M-^) = U-3...(2,-l)^ IZfk' 1 < 30 > 



one finds for the variance of jx 

C{fx) = a 2 a 2 + 2/?V. (31) 
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From fx and the asymmetric errors alone it is not possible to determine the value of the variance, 
since in addition to fx, a and f3 also (x — (x)) is an unknown, i.e. one is faced with three 
constraints for four variables. In the same spirit, however, as fx is an unbiased estimator for fx, 
one can also construct an unbiased estimator C(fx) for C. The solution is given by 

C(il) = \(* + +a-) 2 -l(* + -CT-) 2 . (32) 
One easily verifies that the expectation value of eq.(|32|) reproduces eq. (|3~1~D . The functional 



expression eq. (|32"D can be interpreted as the naive estimate for the variance of fx which is 
reduced by a correction term that vanishes for symmetric errors. An important aspect is that 
C(fx) is only an estimate for the variance of fx, even if the variance a 2 of the measurement 
is known precisely. Only for (3 = 0, i.e. when the connection between measurement x and 
parameter fx is linear one recovers C = C. As shown below, the fact that C(fx) is subject to 
statistical fluctuation complicates matters a lot. 



4.2 Weighted Averages with Asymmetric Errors 

In the previous subsection a prescription was derived to quote an unbiased estimate for mean 
value and variance of a parameter fx. In terms of the expansion eq. (|25|) the results are 

p, = ^ + ay + (3y 2 - pa 2 and C = (aa + 2(3ay) 2 - 2/3 V 4 . (33) 

With eq.(|3~l~l) this can be rewritten as 

C-C = ^(fx-fx) with 7 = 4/5a 2 = 2(a + -a_). (34) 

It follows immediately that the variance of C is given by -y 2 C, i.e. it vanishes for symmetric 
errors. In case of asymmetric errors both C and fx are random variables which scatter around 
their respective expectation values, and according to eq.fl3"4"D both are fully correlated. Although 
individually both variables are unbiased, their ratio fxjC is not. As a consequence also the 
weighted average fx 

fx = — = fx H . (35) 

is no longer an unbiased estimator for the true parameter value fx. Introducing 

di = fxi- fx and d = Q + 7^, (36) 
the expectation value of fx can be written as 

(fx = fx+ { — )=fx+( — ) with Z m = }^— — , 37 

\ E< l /{Ci + lidi) I \Z Q / i c i + idi 

and the bias of the weighted average is given by the expectation value of the ratio Z\/Z$. 

Using a Taylor expansion about the origin, the bias (Zi/Zq) can be expressed through the 
moments of the deviates di. For a general function F of deviates di one has 
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The notation F(0) implies that all derivatives of F are to be taken at di = 0. For the problem 
at hand the first moments are zero, (di) = 0, and different deviates are independent, (d™df) = 
(d™) (df) for k 7^ I. Under these conditions eq.(|38|) simplifies considerably and the leading order 
terms become 

(F\ - I 1 + V {dl) 92 4 V {dl) 93 4 V {dt) ^ 4 V {4) ^ 92 92 I F(0\ (W) 

The evaluation of the moments (df) using eq. fl50p is straightforward. The leading order results, 
expressed as function of the variance Ck and the asymmetry parameter 7/., are 

(d£) = C k , (4) = ^ fc7fe -l 7 3 an d (4) = 3C fc 2 4 3 7 , 2 C fc - (40) 

At this points all ingredients are available which are needed to construct an unbiased weighted 
average also in presence of asymmetric errors. To simplify the organization of terms it is 
convenient to define, in addition to Z m from eq.(^), the following sums: 

For the derivatives of Z m and Ti m , n > 0, one finds 

#m(0) = n! „ n _ m+1 if n > m and ^T /m (0) = ^ n+m — — — . (42) 



dd n k ~" iy "> '"" Cl~ m+l ~ dd n k imw ~ C£ +m (m-1) 

Evidently every derivative picks up one power of 7^ and the expansion around d% = will be 
a power series in 7 j. The error of a truncated series will be quoted as 0(7"), meaning a sum 
where each term is a product of at least n factors 7 j. In general these factors will come from 
different measurements /V 

The actual calculations are straightforward but lengthy. Results of auxiliary calculations are 
collected in the appendix. For the bias on the weighted average one finds 



ZqI Sqi S i 



+ ^+0(7 3 )- (43) 



Equation(|43D shows that in case of asymmetric error the bias is of 0(7), i.e. proportional to 
the asymmetry of the errors. After correction for the dominant effect the residual bias is only 
of 0(7 3 ), because both the third derivative of Zi/Z and the third moments each provide 
one extra power of 7. The higher order terms also contribute at least two factors 7. 

In practical applications a bias correction according to eq.(^) will be done by substituting the 
Ti m for the Si m . This is permitted as long as the expectation value (Tj m ) is equal to Si m to 
the same order in 7 as the terms that are anyhow neglected in the correction terms. As shown 
in the appendix, the additional bias from this substitution for eq. ([43"|) is again of 0(7 3 ). To 
leading order it thus is consistent to use T} m instead of Si m to obtain a bias corrected weighted 
average /2 C 

- , ^11 -^12 , %i , ^11 ^12 A f A .\ 

/i c = /i4 — - 7f ^=fi + — 4 — -=-=// + A (44) 

with p, c — fi = 0(7 3 ). 
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The variance C(// c ) of the bias corrected average \x c is given by 



C(a*c) — (A 2 ) - (A) 2 = (|) -(ft) 2 + 2((ftifc>- (ft) (I 
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+ (t)-(ifc) - 2 «tt>- <*><!» («> 

/?? 2 \ /T 12 \ 2 _ oZ/TnTiaV / Tu\ / T 12 



The individual contributions to this expression can be found in the appendix. Collecting all 
terms one obtains 



C(fe) = ^ + §^-^P + §+0(/). (10 

Doi >J 01 O i O l 



In practical applications one again has to use Ti m instead of Si m , which for the first term in 
eq.(^) introduces an additional bias of 0(7 2 ) which must be compensated by some higher 
order terms. To leading order one has 

-L = i- + ^-^ +0 (/). (47) 

OOl J 01 J 01 J 01 

For all other terms the substitution T} m for Si m is consistent, and the variance C(p, c ) expressed 
through known quantities becomes 



-'Ol J 01 J 01 J 01 



Note that in case of symmetric errors only the first term contributes, recovering the conventional 
estimate for the variance of a weighted average. 

Finally, the expectation value of the x 2 -test-variable shall be examined for a weighted average 
based on measurements with asymmetric errors. Defined in the usual way, \ 2 is given by 

X 2 = E^ l! = ^"f (49) 
k C fc ^o 

where /t*. and are the unbiased estimates for mean value and variance of the individual 
measurements and /2 determined according to eq . (|35|) . Up to second order the determination 
of the expectation value is easily performed, giving the result 

' 2\ _ i i \ ^ S\ x + 45*22 45'ii5'i2 + 75*23 25f 2 n( 4 -> 



X 2 ) = (n - 1) + ^21 - " n ' ~" + ""T ' ^ - ^ + 0( 7 4 ), (50) 

Z O l z >Joi *->01 

where n is the number of measurements contributing to the weighted average. For practical 
applications this result implies, that also in the presence of asymmetric errors the usual \ 2 ~ 
variable can be used as a goodness-of-fit criterion. For a consistent set of measurements the 
value should be not too deviant from the number of degrees of freedom. One should also note 
the somewhat surprising result that a meaningful \ 2 ~ variable is based on the non-bias-corrected 
average p,. 

Before moving on, it is worth while to discuss the explicit expressions when averaging n data 
points which all are drawn from the same parent distribution, i.e. Ck = C and 7^ = 7 for all 
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indices k, as then the inherent structure of the results becomes most evident. One finds the 
following expressions for the expectation values discussed above: 

(fl) = /i- 7 (l-£) + 0(7 3 ) (51) 

(C(fi)) = - + ^-(l--) + 0(7 4 ) (52) 
n n \ nj 

(x 2 ) = (n-l) + ^(n-4 + £) + 0( T 4 ) (53) 

A simple cross check is obtained by setting n — 1 in which case the values (/i) = /i, (C(fl)) = C 
and (x 2 ) = are recovered correctly. One also sees that for 7 = the usual expressions known 
for weighted averages are obtained. The corrections which are needed for non-linear relations 
between measurements and model parameters are an expansion in the non-linearity parameter 
7 and 1/n. 



4.3 A Numerical Example 

Since some of the above results may seem counterintuitive they shall be illustrated by a nu- 
merical example. The model assumes that the relation between a true parameter fi t and the 
expectation value (x) of a measurement x is given by (x) = ^[Jh- Assuming that an experiment 
provides an unbiased measurement x with a Gaussian PDF around the mean, the likelihood 
function is given by 

Pm ~ exp (- {X ~^ )2 ) ■ (54) 

To study the averaging procedure, measurements x were generated in a Monte Carlo simulation 
and treated in the usual way to obtain a Maximum Likelihood estimate fi m i with asymmetric 
errors covering a 68% confidence level interval. The properties of single estimates as well as 
weighted averages will be discussed below. 

For the simulation it was assumed that the true value is fa = 25 and that the variance of 
the measurements is a 2 = 0.5. The characteristics of the problem then are given by (x) = 5, 
a = 10, (5 — 1, 7 = 2 and C = 50.5. A single value x yields a maximum likelihood estimate 
i^mi = x 2 with asymmetric errors a± = 2xa ± a 2 . 

The distribution of the measurements x is shown together with an example for the slightly 
asymmetric likelihood function, taking x = 5, in the top row of fig.([TJ). For this example the 
maximum likelihood parameter estimate becomes fi m i = 25±ggy. The two lower plots of fig.(|l]) 
display the distributions of the maximum likelihood estimates that result from the Gaussian 
distribution in x. The left hand plot shows the distribution of the plain estimate /i m /, the right 
hand plot has the bias corrected distribution of fi. The expectation values are (/i m z) = 25.5 
and (fi) = 25. The plain Maximum Likelihood estimate evidently is biased, consistent with the 
expected bias (cr + — ct_)/2 = a 2 = 0.5. The example nicely illustrates that an asymmetric error 
o + > cr_ really means that the maximum likelihood estimate tends to be larger than the true 
parameter value. To correct for the bias one has to reduce the value by the average asymmetry 
of the errors. 

As an aside, it may be interesting to note that using the mean value (n{x)) of the likelihood 
function /x(x) rather than the maximum as an estimate for the unknown parameter actually 
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Figure 1: Distribution of measurements (upper left), and likelihood function for the 
parameter /i from a measurement x = 5 (upper right). The lower plots show the dis- 
tribution of the uncorrected (left) and the bias corrected (right) Maximum Likelihood 
estimates, [i m i and ft (right), respectively. 



makes things worse. If spill-over to negative values can be ignored, then one finds (/J*(x)) = 
x 2 + 3a 2 , i.e. a bias three times larger than the bias of the Maximum Likelihood estimate \i m \. 
Also this is illustrated in fig.(|l|) by the likelihood function /i(x = 5). 

The Monte-Carlo simulation shows that the variance of the bias corrected maximum likelihood 
estimate C(fi) is consistent with the expectation C = 50.5. It is interesting to compare this 
value to a heuristic estimate for the variance, based on the average of the positive and the 
negative errors. The expectation value of this heuristic variance is slightly larger than C. One 
finds ((a+ + a_) 2 /4) = (Aa 2 x 2 ) = Aa 2 ((x) 2 + a 2 ) = 51. The term (<r+ - oJfjl in eq.® 
corrects this bias. 

Having established unbiased estimates for fi and the variance of that estimate, the next step 
is to study the weighted average of results from independent measurements. As discussed 
before, using the inverse of the variance estimate as the weight, biases the average because of 
the correlation between parameter and variance estimate. The result of using eq. (|35|) for an 
average of 10 values is displayed in the upper left plot of fig.©. The central value of the average 
is about 0.8 standard deviations lower than the true parameter value. The estimated bias is 



shown in the upper right. According to eq . (f43|) , i.e. neglecting terms of 0(7 4 ), one would 



expect a shift of 1.8 units, which is close to what is actually observed. The bias corrected 



12 



weighted averages jHof 10 measurements 




Mean 23.14 
RMS 2.362 





bias estimates b 



5 10 15 20 25 30 35 40 45 p, 
bias corrected averages fi c 



Mean 
RMS 




24.93 
2.368 



10 15 20 25 30 35 40 45 jl c 



Mean 
RMS 



1.782 
0.1371E-01 



0.5 1 1.5 2 2.5 3 3.5 4 4.5 b 
variance estimates for /2 C 



Mean 
RMS 



5.439 
0.4665 



1 



8 10 12 14 16 



18 

c(Ac) 



Figure 2: Weighted average of 10 bias-corrected maximum likelihood estimates fi with 
unbiased variance estimates C(£i). The upper left plot shows that the uncorrected 
average p, is significantly biased. The bias estimate is displayed in the upper right. 
The distribution of the bias corrected averages ]x c is shown in the lower left, and the 
distribution of the variance estimates for the corrected averages on the lower right. 



average is shown in the lower left of fig.(f2|). The lower right finally shows the distribution of the 
variance estimate for the bias corrected weighted average. The distribution is rather narrow 
with a central value close to the actual variance found in the distribution of the averages. That 
the general picture is consistent is also shown in the ^-distribution of the weighted average, 
fig.©. According to eq. (|50"D , including all terms up to 0(7 2 ) one would predict an expectation 
value (x 2 ) ~ 9.25, again very close to the observed value. 

To illustrate further the properties of different schemes to form a weighted average from mea- 
surements with asymmetric errors, fig.(fD shows how the various methods converge when in- 
creasing the number n of measurements. Here n was varied from 1 to 50. The top plot shows 
the result of a simple heuristic method, where the maximum likelihood estimates fi m i are av- 
eraged, using as weights the inverse of (a + + cr_) 2 /4, the middle diagram was obtained by 
using unbiased, though correlated, estimates for mean value and variance, and the lower plot 
shows the result of the proper bias correction scheme. The first two schemes yield very similar 
results, the main difference being a shift corresponding to the bias correction for an individual 
measurement. The heuristic method has a large bias for n — 1, which is corrected in the second 
method. At large n both significantly underestimate the true parameter value. One also has 
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Figure 3: Distribution of x 2 - values for the 
weighted average of 10 bias-corrected max- 
imum likelihood estimates. 



to mention that both provide error estimates which underestimate the actual uncertainty of an 
individual weighted average by roughly 10%. The proper bias correction scheme has both a 
stable expectation value for the average and a reliable error estimate. 



evolution of weighted averages 



heuristic average z 
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Figure 4: Evolution of mean value and 
error estimate as function of the num- 
ber of measurements in the average. 
The top plot shows the behavior of 
a simple heuristic treatment of asym- 
metric errors, the middle one the re- 
sult of using bias corrected estimates 
fx and C{fi) for mean value and vari- 
ance, and the lower plot the results 
after proper bias correction for the av- 
erage and its variance. 
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5 Systematic Uncertainties 

The previous considerations apply to statistical uncertainties. Many measurements are also 
subject to systematic or theoretical uncertainties, which affect the model that is fitted to the 
experimental data. Formally the effect of these uncertainties can be described by an additional 
parameter a entering the model, where the knowledge about a is described by a PDF h(a). 
The function h(a) can be the result of a measurement or a purely Bayesian parameterization 
for the subjective degree of belief in a certain value a. The formal treatment will be unaffected 
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by this. With the parameter a eq.(pID becomes 



1 / (x - s(p,a)) 2 \ 
P{X ^ a) = ^Ta 6XP 2^ • (55) 



The simplest and probably still very common approach to deal with the situation is to vary 
a according to h(a) around a nominal value ao over a range covering a 68% confidence level 
interval, and quote the resulting variation in p as a systematic error. The discussion presented 
here will deal with only this approach. As before, from the quoted errors the mean value and 
variance of the parameter estimate shall be inferred, since that is the information which is 
needed to form a weighted average. 

For a fixed parameter a the maximum likelihood estimate p m i for p according to eq.flSSD is 
given by p m i = s _1 (x, a). The resulting PDF p(p m i) for p m i depends on a) and the PDF 

h(a) of a. 

p(p m l) = J 'dah(a) 5 (p m i - s~ l (x,afj . (56) 
From eq. (|56|) one finds for the moments of p m \ 

Vml) = f dah(a) (s^(x, a)) k . (57) 



To continue, it will be assumed that h(a) is a Gaussian with mean value ao and standard 
deviation r, and that s~ 1 (x, a) has an expansion around ao 

s (x, a) = po(x) + a(x)(a — a ) + /?(x)(a — a ) 2 (58) 

where the higher order terms are negligible. To simplify the notation the x-dependence of the 
coefficient will be dropped. With 

%)= 1 expp^l (59) 



one finds for the first and second moment of p m i 

(Pmi) = Pa + (3t 2 and (p 2 ml ) = p\ + (a 2 + 2p (3)r 2 + 3(3 2 r\ (60) 



From eq.(|60|) one finds the variance of p m i 

C{p ml ) = (p 2 ml ) - (p ml ) 2 = a 2 r 2 + 2/3V. (61) 

For a = ao one has the nominal parameter estimate «o- The in general asymmetric errors are 
given by 

t + = ar + j3r 2 and r_ = ar - j3r 2 . (62) 

Combining the information from eq.(^) and eq.(|B"2"D the unbiased parameter estimate and its 
variance can be written as 

(Pmi) = Po + ^(r + - r_) and C(p m i) = ^(r + + r_) 2 + ^(r+ - r_) 2 . (63) 

These expressions eq.(|63"D show how to deal with theoretical or systematic uncertainties when 
averaging independent measurements. In contrast to the statistical errors discussed before, now 
asymmetric errors translate into a shift of the parameter as one would naively expect. Also the 
variance is strictly positive. Although the derivation is formally very similar to the one of the 
statistical errors, here the PDF of the function of the parameter a is known a priori, which is 
not the case for the measurement x. Note that the expressions eq. (|63"D have been derived for 
the assumption that systematic uncertainties are Gaussian errors. 
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6 Combination of Statistical and Systematic errors 



Except in the context of Bayesian statistics the combination of statistical and systematic errors 
is not a well defined concept. If one, however, adopts the attitude to model the influence of 
the two kinds of uncertainties by means of a Monte Carlo method, where the measurement 
x as well as the theory-parameter a is varied according to their respective probability density 
functions, then the total bias and variance for the parameter estimate ft is given by the sum of 
the individual biases and the individual variances, i.e. the errors can be combined in quadrature: 

Ck = Cf at + Cf st . (64) 

In the Bayesian sense this is the correct way of combining different sources of uncertainty. From 
a non-Bayesian point of view the ansatz eq. (|64]) at least has the qualitatively desired properties, 
that measurements with large uncertainty in either of the two components will get a low weight 
in any average. 

When using the combined variances in a weighted average, one again needs to perform a bias 
correction based on the asymmetry parameters 7. As defined in eq.(^4|), only the difference 
of the statistical errors contributes to 7, since only asymmetries in the statistical uncertainties 
lead to a correlation between variance and mean value. No such correlation exists for systematic 
errors. Otherwise the definitions eq.(STj) apply as before. 



7 Summary and Conclusions 

The weighted average with weights based only on the variance of the individual measurements 
is a robust estimator for a common mean. As a consequence of the central limit theorem, its 
PDF can be assumed to be a Gaussian even if the PDFs of the inputs are not known in detail. 
Although it would be preferable from a theoretical point of view to take correlations into account 
when forming the average, this can only be recommended if those are well understood. A simple 
effective way is proposed to include the impact of correlations between the measurements for 
the error estimate, using the \ 2 °f the average as an indicator variable. 

In a second step the problem of asymmetric 68% confidence level error intervals was addressed. 
The argument was based on the assumption of a measurement with a Gaussian PDF and the 
use of the Maximum Likelihood Method to infer the value of a physical parameter and its error. 
Explicit expressions were derived to convert the given information into an unbiased estimate for 
the value of the parameter and its variance. In addition, it has been shown that when dealing 
with asymmetric errors, correction terms for the weighted average and its variance have to be 
taken into account. In case of symmetric errors the usual expressions for the weighted average 
are recovered. The scaling prescription to account for correlations between measurements can 
be applied in the same way for measurements with symmetric and asymmetric errors. 

Finally, systematic errors were studied, assuming that they enter as model parameters with a 
known PDF that has to be interpreted in the Bayesian sense. As such they behave differently 
from the statistical errors of a measurement. Also here explicit expressions are given how to 
infer mean value and variance and how to combine them with statistical uncertainties. 
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Appendix: Auxiliary Expectation Values 

Here the complete set of expectation values is given which are needed to derive the results of 
section |4~2] . In most cases a precision of at least 0(7 2 ) for the expectation values is obtained 
already when truncating the expansion eq. (|39|) after the second term. 



<^> = n0)+£^!^(0) + ... (65) 
, z. ua k 



Elementary calculus with the derivatives and definitions given in eq. fl3"7j) , eq.fl41|) and eq. (142|) 
then yields the following results: 



-7T- ~ 7? + S + °(7 4 ) (66) 



1 \ _[ S22 S23 n , _ 1 

S01 S 01 S 01 
Tn\ Sn : S32 S11S22 + 5*33 : S11S23 n , 5 



T-l / - Q + Q ^ + ^3— + °VT) ( 67 ) 

J 0l/ >->01 >->01 ^01 D 01 

/ T n \ Sf x 2S11S32 + 5*43 25^ 5*22 + 4S'ii5'33 SS^S^ „, 6 > , N 

(^T/ = OF + C2 ^ + C4 + °(7 ) (68) 

\ J 0l/ °01 °01 °01 D 01 

(7?) = ^ + ^ + ^r^+0( 7 5 ) (69) 

\ J 01 / °01 °01 °01 °01 

/ ^12 \ ^12 1 6S , i 2 5 , 33 + 4S45 4S ,2 2 S , 2 2 + 16S , i 2 5 , 34 lOS^S^ „, 6N , 

(^r> = + ^ ^5 + — ^ — + 0{Y) (70) 

\ J 0l/ °01 °01 '-'Ol °01 

Ti\T\2 \ 5*ii5'i2 S12S32 + 35 l n5 l 33 + 26*44 ^ 7 ,n 

, J 01 / °01 °01 

35*11 S'i25'22 + 3S'i2>S'33 + 6S'iiS'34 , 6S'ii5'i2S'23 , 



C4 1 C5 

°01 °01 



+ — + o( 7 b ) 



^0/ 5"oi S m 



+ 0( 7 3 ) (72) 



TnZi \ 5i 2 i + S22 + 25^2 + 0(74) (73) 



= - S " S 'L +2S23 + ^ + 0(y) (74) 
v i oi z, o/ <jqi D oi 
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The quantities where one has to go to fourth order in eq.([39|) in order to collect all terms up 
to 0( 7 2 ) are (Z^/Zq), (Zf/Z ) and (Z 2 ). For F = Z\jZ^ the partial derivatives up to fourth 
order are: 

J_ Wr* 1 2 7fc7; 1 / 2 7 2 4 7fcT; 4 TfcT; 27? \ 

2!2!«9d 2 cM 2 W " S&CgC? SiyClCf^ ClCf + ClCf + ClCf) ^ } 

, 1 f 3 7 2 12 7fc7 , 3 7 2 ^ 



c4 I riir~i2 r~i3r^3 ry2ri4 I 
D 01 K^k^l ^k^l ^k^l J 

With the moments given in eq. (|^0|) and using the identity 

2 £ a fc 6, = ( £ aj M - fe ( 79 ) 

fc<Z \fc / \ l ) \ k ) 



the expectation value (Zf/Z^) finally becomes 



Z l \ 1 5 n + 35 2 2 4^n^i2 + 6^23 65i 2 h 

■^2 / = c~" H Q2 C3 h "^4" + <A7 )• 

Z/q i Ooi Oqi >J i Oqi 



Similarly one finds for the other expectation values 



ZA = 1 + ^i±i^ - ASll %t 7823 + 2 ~¥+ Q(7 4 ) (81) 
(^2) = ^ + ^2i + 0( 7 4 ) (82) 



where n is the number of measurements contributing to the average. 



18 



